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ABSTRACT 

We discuss the thermodynamics of gluons in the background of static quark 
sources. In order to do so we formulate the quenched limit of QCD at non-zero 
baryon number. A first numerical analysis of this system shows that it undergoes 
a smooth deconfining transition. We find evidence for a region of coexisting phases 
that becomes broader with increasing baryon number density. Although the action 
is in our formulation explicitly Z(3) symmetric the Polyakov loop expectation value 
becomes non-zero already in the low temperature phase. It indicates that the heavy 
quark potential stays finite at large distances, i.e. the string between static quarks 
breaks at non-zero baryon number density already in the hadronic phase. 
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1 Introduction 



Numerical studies of the SU(3) gauge theory, i.e. the heavy quark mass (quenched) 
limit of QCD, are extremely helpful for the understanding of the phase structure 
of QCD at non-zero temperature. The deconfming phase transition as well as basic 
properties of the low and high temperature phases can be studied in this approxi- 
mation, although the phase transition in the light quark mass (chirat) limit differs 
not only on a quantitative level but also qualitatively; in the chiral limit the order 
of the transition is controlled by the chiral symmetry of the fermionic part of the 
QCD Lagrangian whereas in the quenched limit the Z(3) center symmetry dictates 
the order of the transition. Nonetheless, the fundamental features of deconfmement 
which are reflected, for instance, in the release of many new degrees of freedom at 
T c , are similar in both limits. This leads to a quite similar temperature dependence 
of bulk thermodynamic observables like the pressure and energy density 

So far the investigations of QCD thermodynamics were essentially limited to the 
case of vanishing baryon number. In the case of non-zero baryon number, usually 
realized in thermodynamic calculations through the introduction of a non-zero chem- 
ical potential (//) 0, ||, little is known from lattice calculations about the behaviour 
of thermodynamic observables, the QCD phase diagram and the properties of the 
different phases. The reason for this is well-known. The probabilistic interpretation 
of the path integral representation of the QCD partition function breaks down as 
soon as one introduces a non-zero chemical potential. Moreover, even the naive 
quenched limit at non-zero chemical potential, i.e. the ordinary SU(3) gauge theory 
at fi = 0, turned out to be pathological when fermionic observables with /i ^ 
are analyzed M. In fact, it is understood that this naive limit is not the correct 
quenched limit of finite density QCD; it is the zero flavour limit of a theory with 
equal number of fermion flavours carrying baryon number B and —B, respectively. 

As the relativistic chemical potential \x also contains a contribution proportional 
to the rest mass of the quarks, the static limit of QCD at non-zero chemical potential 
requires to take a double limit in which the quark mass as well as the chemical 
potential is taken to infinity while an appropriate ratio is kept fixed. This limit 
has been formulated in || f7|. It seems that in this case the first order deconfming 
phase transition of the SU(3) gauge theory turns into a crossover for arbitrarily 
small, non-zero values of the chemical potential [0]. However, although at non-zero 
chemical potential the thermal phase transition may get lost as a true singularity 
of the partition function of quenched QCD|, we still expect the quenched theory to 

a This may even be the case under more general circumstances. If the phase transition at finite 
density is more like a percolation transition thermal observables may be non-singular although 
singularities due to percolation of domains with high energy or particle density will still exist 
(Kertez line §) [§. 



2 



resemble the physics of deconfmement - gluon thermodynamics in the background 
of a non-zero number of static quark sources is expected to lead to deconfmement 
at high temperature. It is then interesting to analyze the nature of this transition 
and study in how far static quark sources in a gluonic heat bath influence the 
deconfmement of gluons and, for instance, the heavy quark potential. We will 
discuss here a framework for the analysis of these questions - quenched QCD at 
fixed baryon number - and will present first numerical simulations which address 
these questions. 

Rather than introducing a non-vanishing chemical potential, i.e. formulate QCD 
at non-vanishing baryon number density in the grand canonical ensemble, one may 
go over to a canonical formulation of the thermodynamics and fix directly the baryon 
number jlOfl . This is achieved by introducing an imaginary chemical potential [K], 



TT|| in the grand canonical partition function. Performing subsequently a Fourier 
integration allows to project onto the canonical partition function for a given sector 
of fixed baryon number |10| |. In the heavy quark mass limit one is then left with 



a well-defined quenched theory at fixed baryon number - gluon thermodynamics in 
the background of a non-zero number of static quark sources, suitably arranged to 
obey Fermi statistics. 

Despite the fact that the grand canonical formulation of finite density QCD 
leads to severe numerical problems ]12| the alternative canonical approach at non- 



zero baryon number so far did find only little attention^. To some extent this 
may be justified. Introducing a constraint on the total number of fermions (quarks 
minus anti-quarks) in general leads to rather complicated non-local constraints for 
the gauge field sector. However, in view of the problems that arise in the standard 
non-zero chemical potential formulation and that grow exponentially with the size of 
the volume one may seriously want to analyze also the canonical approach in more 
detail. We will show here that the formulation of QCD at non-zero baryon number 
does have a well defined, non-trivial heavy quark mass (quenched) limit which can 
be analyzed numerically at least for moderate baryon number densities. 

In the grand canonical formulation the Fermion determinant becomes complex 
for any non-zero chemical potential leading to all the conceptual and algorithmic 
problems. In the canonical approach, on the other hand, the determinant stays real. 
However, as should be obvious, the problems related to a non-positive integration 
measure are not solved that easily. The Fourier transformation which is needed 
in the canonical formulation to project onto the sector with fixed baryon number 
reintroduces negative contributions to the partition function and one thus again faces 
algorithmic problems. Performing the Fourier integration numerically thus seems 



b It seems that in addition to the mean-field analysis performed in Ref. |l0| ] in the heavy quark 
mass limit so far the canonical formulation only has been used in a numerical study at strong 
coupling (/3 = 6/g 2 — 0) using the Monomer-Dimer-Polymer representation of QCD (l3[] . The 
potential power of this approach has, however, also been stressed recently in Ref. p~4|| . 
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to be ruled out; one can, however, do it explicitly. This leads to a complicated 
expression in terms of products of quark propagators which may be viewed as the 
boundary conditions for the gauge fields needed to project onto a sector with fixed 
baryon number. We will use this as a starting point for our analysis of the heavy 
quark mass, quenched, limit. 

In the next section we will discuss the canonical formulation of QCD at finite 
baryon number and, in particular, introduce the quenched limit. In section 3 we will 
present some results of a first numerical analysis of this quenched theory. Section 
4 contains our conclusions. In an appendix we give explicit representations for the 
canonical partition function for some values of the baryon number. 



2 Lattice formulation of QCD with non-zero 
baryon number 

The general framework for going from a grand canonical formulation of QCD at 
non-zero baryon number density in terms of a non-vanishing chemical potential 
0, |3|] to a canonical formulation with a fixed baryon number has been outlined in 
[10|| . As mentioned before the main difficulty of this approach arises from the need 



to perform a Fourier transformation to eliminate the chemical potential in favour of 
a fixed baryon number. We will show in the following that this integration can be 
performed explicitly. The projection onto a given sector of fixed baryon number is 
then contained in a rather complicated sum over products of quark propagators. This 
representation probably is too complicated to be useful for numerical calculations 
with arbitrary, e.g. the physically interesting, light quark masses. However, it may 
provide a new starting point for physically relevant approximations to the complete 
problem. In particular, we will discuss in the next section the quenched limit at fixed 
baryon number, i.e. gluon thermodynamics in the background of static quarks. 

Starting from the standard formulation of QCD at non-zero chemical potential 
P] and the corresponding formulation at non-zero baryon number [JO] we will rewrite 
the fermion sector of the lattice action in a somewhat more transparent form which 
makes clear that the non-vanishing chemical potential can be viewed as a modifi- 
cation of the temporal boundary conditions for the fermion fields. To be specific 
we will use the Wilson formulation of the fermion action. Although we will eventu- 
ally restrict our discussion to the static limit, which also can be obtained directly 
from the heavy quark formulation for Wilson fermions derived in analogy to the 
approach given in Ref. for staggered fermions, we will start here by formulating 
the canonical approach for arbitrary quark masses. 

The action for Wilson fermions at non-zero chemical potential is given by 
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x \ j=l 

-«[e^4(l - 74)^,4^+4 + e" Ma ^ + 4(l + 74)^,4^]) • (2-1) 

Here k is the hopping parameter which controls the value of the quark mass, 
x = (x, X4) denotes the sites on a lattice of size N% x iV T . The Grassmann fields 
obey anti-periodic boundary conditions in the time direction (fourth direction), 
i.e. ip(x,N T +i) = -ip(3,i) and ^(x,n t +i) = -$(x,i)- 

We may shift the dependence on the chemical potential to the last time slice, 
which connects the hyperplanes with X4 = N T and £4 = 1, by transforming the 
fermion fieldsf] 

^(x,x 4) = ^ ax ^ ( x,x 4) , ^ >4) = e-^ {2>4) • (2.2) 

This eliminates the dependence on the chemical potential on all time slices but the 
last one, which now reads 

S^Qi/T) = /c£[e^ (W (l-74)^ , (2.3) 

x 

with x = (x,N T ). We also have used the definition of the temperature 1/T = aN r 
in writing \ijT — fiaN T and furthermore explicitly took care of the anti-periodic 
boundary conditions for the fermions. The remaining part of the fermion action, 
which now is independent of the chemical potential, may be written as 

S F = S F (0)-S^(0) . (2.4) 

This representation explicitly shows that the formulation of thermodynamics 
with non-zero chemical potential can be viewed as a generalization of the non-zero 
temperature case, which is realized through anti-periodic boundary conditions in 
the fermion sector. They are now generalized to ip(g,N T +i) — — exp (/x/T) V^i) and 
fcjv T +i) = -exp(-/i/T) ^x). 

c Note that the Jacobian of this transformation equals one. 
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So far we have only re-organized the various terms in the fermion sector of the 
QCD partition function where we have used the standard Wilson formulation. For 
the gluon sector we do not have to specify at this point the explicit form of the action, 
Sq. The partition function in a volume V = (N a a) 3 at temperature T = l/N T a and 
non-zero chemical potential /ia then reads, 

Z(/i/T,T,V) = |n d ^n d ^d^e- s " r ^e- 5G -^ . (2.5) 



We want to use this form as a starting point to go over to a formulation at non-zero 
baryon number rather than at non-zero chemical potential. This can be achieved by 
introducing a complex chemical potential and performing a Fourier transformation, 

1 r 2lT 

Z(B, T, V) = — / d0 e' iB * Z{i<f>, T, V) , (2.6) 
Ltx Jo 

where B denotes the quark number, i.e. the baryon number equals B/3. The Fourier 
transformation only operates on the factor e~ f which only involves links pointing 
in the 4th direction on the last time slice of the lattice. Making use of the Grassmann 
properties of the fermion fields this contribution can be written as 



e-^ T ™ = n 1 - ™«^ T ^m >< 

(x,a,b,a,f3J) V 7 

n ( i - ™-w£pp aA "^ii^ , (2.7) 



(x,a,b,a,f3,f) 

where the product runs over all possible combinations of indices with x taking values 
on the three dimensional (spatial) lattice of size N%. We also have introduced the 
notation Ug = T-U^g^)^ and iftg = F + l/L N ^ 4 with T± — (1 ± 74). Note that the 
fields U carry spinor and color indices which we denote by Greek and Latin letters, 
respectively. We will combine these to a single index denoted by, e.g. A = (a, a) 
with a = 1,...,4 and a — 1, 2, 3. In addition we also have allowed for different 
fermion flavours, / = 1, ...,n/ but ignored the possibility of having different quark 
masses, i.e. different hopping parameters k for the various flavours. 

We may expand the product appearing in Eq. [2.7| and write it as a series in 
terms of the complex fugacity, z = exp(z0). The Fourier transformation in Eq. |2l] 
will receive a non-zero contribution only from the term proportional to z B . The 
coefficient of z B will receive contributions from n terms proportional to z and n 
terms proportional to z* where B = n — n. Each such contribution is proportional 
to K n+n . This is the basis for a systematic hopping parameter expansion for the 
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boundary term. Clearly we need at least B terms proportional to z. The leading 
contribution thus arises from the n = sector. It can be summarized as 

z B f B = {-zK) B £ U^.i4: n c n ,;; J u , (2.8) 

X,C,T>,Fi=l 

where X, C, V, F are B- dimensional vectors, i.e. X = (xi, xb), F = (/i, fs) 
and so on. Of course, all elements of the set of indices {(Cj, fi, Xi)}f =1 as well as 
{(T>i, fi, Xi)}f =1 have to be different to give a non- vanishing contribution to the sum 
in Eq. |2.8| . The Fourier integral in Eq. |2]6| can then be performed explicitly and we 
obtain for the partition function at fixed baryon (or quark) number 

Z(B, T, V) = f I] dU x , v J] d$ x di/> x f B e- SG ~^ . (2.9) 

J X,U X 

To leading order in the hopping parameter the projection onto a sector with 
fixed quark number B is encoded in the function fs which is a sum over products 
of quark propagators between the time slices at x 4 = 1 and x 4 = N T . In higher 
orders one will, of course, also get contributions from quarks propagating backward 
in Euclidean time. In fact, if we think in terms of a hopping parameter expansion 
(heavy quark mass limit) for the entire fermion determinant, the function J'b is all 
we need to generate the leading contribution, which finally will be 0(k BNt ). Higher 
order contributions will result from ^-independent terms coming from an expansion 
of exp(— Sf) as well as from additional factors in the expansion of Eq. |2.7| which then 
have to contain an equal number of additional backward and forward propagating 
terms. 

Let us look in more detail at the leading contribution arising from For this 
it is convenient to evaluate Jb, which of course is a gauge invariant function, in 
a specific gauge. Let us perform a gauge transformation such that all the links 
pointing in the time direction on the last time slice are equal to unity. This gives 

fB = (-2*) B £ n^/,i)<i'i) • ( 2 - 10 ) 

X,A,F i=l 

Like in Eq. |2.8| the vectors X, A, F are of length B. However, now the spinor 
indices «j which are part of A% only take on the values a« = 1, 2 because only the 
two components of T_ are non- zero. This also gives rise to the factors of 2 in front of 
k. When evaluating the Grassmann integrals each of the ip terms can be contracted 
with all those if) terms which carry the same flavour index. Each pair gives rise to 
a matrix element of the inverse of Q, the fermion matrix corresponding to Sf- The 
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different pairings give rise to the Matthews-Salam determinant [jnj. We thus will 
get the product of rif determinants, each of dimension d\ such that Ya=i di = B, 



f B = (2K) B l[detMi[x,A) , (2.11) 

X,A,F 1=1 

where the matrix Aii gives the contributions for the 1-th flavour and the matrix 
elements are the corresponding quark propagators, 

M-]' 3 = 0(4,i),^),((^,a^),A) • (2-12) 

Each matrix element of Aii is 0(K^ NT ~ 1+ ^ Xi ~ x ^). In the limit of small k- values, 
i.e. in the heavy quark mass limit, only matrix elements of Ai with |x, — Xj\ = 
will thus survive. In this case the elements of Q7^>. u a a tv t ) a) are J us ^ P rocm cts 
of terms T_U(£ t) k),4 with k — 1, .., A^ T — 1. As T_ is a diagonal matrix in spinor space 
the indices and aj have to be identical. The spinor part thus gives rise to an 
overall factor 2 JVt_1 for each i, i.e. we obtain B such factors. The multiplication of 
the SU(3) matrices yields an element of the ordinary, complex valued (!) Polyakov 
loop (U = 1 on the last time slice !) which we denote by L^' aj . Finally, the sum 
over different colour indices appearing in Eq. [2.11| leads to contributions involving 
only traces over powers of the Polyakov loop, 



Lx= II %,*4) • (2-13) 

£4 = 1 

As the (colour, spinor) label Ai can take on six different values the determinant 
is non-zero only if at most six quarks of a given flavour occupy a given site Xi. 
At most three quarks can have the same spinor component. There are thus six 
possible contributions, D n , of a given site to the determinant depending on the 
quark occupation number, n, of this site. These can be expressed in terms of three 
functions Mj which correspond to the number of quarks (i) with identical spinor 
components. These contributions are 



D 1 = 2Mi 

D 2 = 2{M\ + M 2 ) 

D 3 = 2 i 3.1 / ; .\/, + M 3 ) 

L> 4 = 2(4MxM 3 + 3Mf ) 

D 5 = 20M 2 M 3 

As = 2OM3 2 , (2.14) 



8 



with 



Mi = TrLg i , M 2 = (TiLgJ 2 — TrL§. , M 3 = 6 . (2.15) 

We note that this representation is quite similar to that of the heavy quark mass 
limit for QCD with staggered fermions at non-zero chemical potential ]/]]. 

We now may write fs as 

B' rif 

f B = (2^11111^(4) , (2.16) 

X,Fk=l 1=1 

where B' < B is the number of distinct sites x k appearing in X and n k j is the 
occupation number for these distinct sites x k with quarks of flavour I. They obey 
the constraint J2 k i n k,i — B with < n k ,i < 6. The total number of quarks at the 
site Xk is nk = J2i n k,i- The sum over all different distributions of the flavour number 
can be performed locally for a given distribution of sites X. Let us define 



n k n k 



Tit rif 



Dn k (x k )= 8( n k-J2n kj i)]jD nkil (x k ) , (2.17) 

"■k,l=0 nk,n f =0 1=1 1=1 



with n,k = min(nfc, 6). We then can write Eq. |2.16| as 



f B = (2 K ) BN ^f[D nk (x k ) ■ (2-18) 

X k=l 



One of the difficulties with using a representation like Eq. 2.16 or Eq. 2.18 in 
a numerical calculation is that there appears a -B-fold sum over the volume, V, 
with constraints on the occupation number for a given site, i.e. the computational 
effort would be 0(V B ). Through a cluster decomposition we can, however, reduce 
this to the evaluation of certain moments of Dj as well as their average on a given 
configuration. This is described in the Appendix. 



We finally obtain 



ft 



6nt 



(2k)^£5 fi-EE haai 

{g a } V 1=1 {a} 

1 r 6 / 6n / r 6n / 1 

n fr((-i) (E - ai - 1] {^ - 1)1 II rrA ai 



{«} 



1=1 



1=1 



ail 



fja 



(2.19) 



9 



where [...] is the mean value defined in ( |A.6| ) and the product runs over all possible 
sets of vectors a — (a\, a 6nf ) which are constrained by 



6nf 

Y,lai<B . (2.20) 
1=1 

The sum over g a is defined in ( |A.11| ). 

Now we have achieved that only simple sums over products of suitably chosen 
terms D n have to be evaluated to obtain the contribution fs to the partition func- 
tion which results from a given number of static quarks. The number of possible 
partitions contributing to Eq. |2.19| grows, however, rapidly with B. After all it is 
nothing else but the result of an explicit evaluation of the fermion determinant in a 
fixed baryon number sector. As fs contains products of mean values calculated on 
a given configuration, it is a non-local quantity which in a Monte-Carlo simulation 
has to be evaluated for each update of a temporal link. In practice, calculations 
thus will be limited to small values of B. Some explicit representations of fs for 
small values of B are given in the appendix. In the next section we will use this rep- 
resentation of the quenched limit of the QCD partition function with fixed baryon 
number. We will perform first exploratory Monte-Carlo simulations to analyze the 
phase structure of gluon thermodynamics in the presence of static quarks. 



3 Simulation of quenched QCD with non-zero 
baryon number 

For any fixed value of the baryon number we can write the quenched partition 
function as 

Z(B,T,V) = J l[dU XjV f B e~ s ° , (3.1) 



where the constraint on the baryon number is encoded in the function fs given in 
Eq. ^.lflj . In particular, we note that the global Z(3) symmetry of the QCD partition 



function at non-zero baryon number is preserved also in the quenched limit, i.e. the 
function fs is invariant under global Z(3) transformations if B is a multiple of 3. As 
the gluonic action Sq also shares this property the partition function, Z(B,T,V), 
is non-zero only if B is a multiple of 3. 

We also note that fs is still a complex function. However, upon integration over 
the gauge fields the imaginary contribution vanishes; the partition function is real, of 
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course. Actual calculations thus can be performed with Kefs- The crucial question 
for a simulation with this partition function is whether the remaining sign changes 
of the real part of fs are seldom enough so that a Monte-Carlo simulation can be 
performed with the absolute value of Re/s and the overall sign can be included in 



the calculation of averages [J| 



We will perform simulations for the one flavour case (rif = 1) using the partition 
function 

Z n (B,T,V) = / J] dU XiU \Ref B \ e~ SG . (3.2) 



Expectation values of an observable O calculated with the Boltzmann weights de- 
fined by the partition function Z(B,T,V) will therefore be calculated according 
to 

= (0-sgn(Re/ B )>„ (3 3) 

(sgn(Re/ B ))|| 

Our simulations are performed on lattices of size 8 3 x 2 and 10 3 x 2 using the standard 
Wilson gauge action |T7|]. For the link updates we use a Metropolis algorithm. 
For each link update the change in the function Kefs is calculated and a possible 
change in sign is monitored. In all the cases we have studied so far we find that 
(sgn(Re/s))|| is large and can be well determined. In fact, for large values of the 
temperature sgn(Re/e) is almost always positive. This is evident from Figure p] 
which shows the average sign as a function of the coupling j3. The expectation value 
of sgn(Re/e) mainly depends on the spatial volume N a but varies little with B at 
fixed N a . We also find that the values of observables like the average action or the 
Polyakov loop do not depend much on the sign of Refs so that the errors obtained 
for these observables from a jackknife analysis are substantially smaller than those 
shown in Figure [I]. 

A numerical analysis of the thermodynamics at fixed baryon number, B, can 
closely follow the standard approach at B = 0, i.e. in a pure SU (3) gauge theory ||I$f . 
We may analyze the temperature dependence of bulk thermodynamics, the Polyakov 
loop expectation value and other observables for a gluon gas in the background of 
static quarks. We started a first exploratory analysis of this system by performing 
a numerical simulation on lattices with temporal extent N T = 2. The simulations 
have been carried out in the vicinity of the critical coupling for the deconfinement 
transition at B — 0, i.e. for gauge couplings f3 ~ 5.0 which are still in the strong 
coupling regime. Calculations with fixed B are performed on lattices of size iV 3 x N T 
and the temperature is varied, as usual, by changing the coupling (3 = 6/g 2 . The 
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Figure 1: Expectation value of the sign of Re/s, (sgn(Re/s))||, for B 
and lattices of size iV 3 x 2 with N„ = 8 and 10. 
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dimensionless parameter kept fixed in the simulation thus is the baryon number 
density in units of the temperature cubed, 



T 3 3 \N a 



(3.4) 



The baryon number density in physical units thus is 



n B 



3 \N a 



T 



200 MeV 



fm 



(3.5) 



For orientation we note that close to T c , which for the SU(3) gauge theory is known 
to be about 270 MeV, a simulation on an 8 3 x 2 lattice with 5 = 12 corresponds to 
tib — 0.15/fm 3 , i.e. approximately nuclear matter density. 

For vanishing baryon number the Polyakov loop expectation value or more pre- 
cisely the expectation value of its normalized absolute value calculated on a finite 
lattice, 



W>v = <I™Em*I> 



3iV 3 _ 

o x 



(3.6) 



is an order parameter for the deconfinement transition in the infinite volume limit, 



(L) = lim (\L\) V 

N a —>co 



(3.7) 
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As the phase transition is first order in 577(3), (L) changes discontinuously at T c . 
Besides being related to the spontaneous breaking of the global Z(3) symmetry of 
the pure gauge action the behaviour of the Polyakov loop also reflects the large 
distance behaviour of the heavy quark potential 

e - V {x-y,T)/T = (TrL^TrLt) — > 9|(L)| 2 . (3.8) 

\x — y\ — > oo 

The vanishing of (L) indicates that the heavy quark potential is confining at large 
distances, while a finite value of (L) shows that the potential stays finite for infinite 
separation of the quark anti-quark pair. In QCD with dynamical light quarks the 
Polyakov loop is no longer an order parameter. The heavy quark potential stays 
finite at large distances even in the confined phase because the static quark anti- 
quark pair can be screened through the creation of a light quark anti-quark pair 
from the vacuum (string breaking). For light enough quarks this is indeed observed 
in numerical simulations at vanishing baryon number density |18 . 



At non-zero baryon number density we expect to find a similar behaviour of the 
heavy quark potential even in the heavy quark mass limit because the quarks needed 
to break the string need not be created through thermal (or vacuum) fluctuations. 
The static quark anti-quark sources used to probe the heavy quark potential can 
recombine with the already present static quarks and will lead to string breaking 
even in the low temperature hadronic phase. We thus expect that the Polyakov 
loop expectation value will not be an order parameter, although the integrand of 
the partition function, /sexp(— Sg), is Z(3) symmetric, i.e. we expect that 

(L) > , for all n B > and all (3 > . (3.9) 

This is indeed evident from the results obtained for the Polyakov loop expectation 
values from our simulations with B = 6 and 12 on 8 3 x 2 and 10 3 x 2 lattices which 
are shown in Figure [| We thus find first (indirect) evidence for the modification 
of the long distance part of the heavy quark potential in nuclear matter. This will 
be analyzed in more detail in the future. We also note that there is no significant 
volume dependence at fixed ub^\ This also shows that in the thermodynamic limit 
the physical observables will indeed only depend on the density rather than the 
baryon number itself. 

For B = there is a clear signal for a first order phase transition which leads 
to a discontinuity in (L). For all B > we clearly observe a transition from a 



Calculations performed on a 8 3 x 2 lattice with B = 6 and on a 10 3 x 2 lattice with B = 12 are 
performed at nearly the same baryon number density, i.e. ns/T 3 = 0.03125 and ns/T 3 = 0.032, 
respectively. 
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Figure 2: Polyakov loop expectation value for different values of B and lattices with 
spatial extent N a = 8 and 10. 



low temperature phase with small Polyakov loop expectation value to the high tem- 
perature regime characterized by a large Polyakov loop expectation value, which is 
similar to that of the B = case. The transition occurs in a temperature interval 
that broadens with increasing baryon number density. There is no indication for a 
discontinuous transition. In fact, this is not to be expected in a canonical calcu- 
lation, even if the transition is first order. By changing the gauge coupling /3 we 
vary the lattice cut-off and through this also the baryon number density continu- 
ously. At fixed non-zero baryon number we therefore follow a simulation path that 
traverses continuously through a region of two coexisting phases. This situation is 
schematically illustrated in Figure |3|. 

The question now is whether the transition region really is a region of coexisting 
phases. In this case the values of thermodynamic observables result as the super- 
position of contributions from two different phases appropriately weighted by the 
fraction each phase contributes in the coexistence regime. In the infinite volume 
limit this could result in a discontinuous change of the slope of thermodynamic ob- 
servables when entering and leaving the coexistence region (for an illustration see 
e.g. Figure [3|b). 

To gain further insight into the structure of this regime we also analyzed various 
susceptibilities. In Figure |] we show the conventional Polyakov loop susceptibility, 



XL = iV ff 3 ((|L| 2 )-(|L|) 2 ) 
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(3.10) 
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Figure 3: Schematic plot of the QCD phase diagram (a) in the temperature-baryon 
number density plane for the case of first order transitions in the entire plane. For 
B > and T h < T < T q the system stays in a region of two coexisting phases. For 
B — the transition occurs at a unique temperature T c . In (a) we also show the 
paths followed when varying the coupling (3 in a Monte-Carlo simulation with fixed 
B, N a and N T . In (b) the expected behaviour of the Polyakov loop expectation 
value along these paths of non-zero B as well as for B = is shown. 



as well as the derivative of (\L\) with respect to (3, 
d(\L\) 



Xf3 



dp 



(\L\.S G )-(\L\)(S G ) 



(3.H) 



Both response functions reflect the existence of a transition region that becomes 
broader with increasing Compared to the behaviour at B = they also change 
continuously in this region. Such a behaviour might as well just correspond to a 
smooth crossover to the high temperature regime; a conclusion also drawn from 
the heavy quark simulations with non-zero chemical potential J7|. To establish the 
existence of a region of coexisting phases with certainty will thus require a further, 
detailed analysis of finite size effects. 

The width of the transition region gives some indication for the shift of the 
critical temperature when going from B = to a baryon number density which 
roughly corresponds to nuclear matter density. The gauge coupling at which the 
Polyakov loop expectation value starts rising rapidly is shifted from (3 C = 5.09 at 
ub/T 3 = to (3 C ~ 4.95 at hb/T 3 = 0.032. A rough estimate based on the non- 
perturbative /3-function for the SU(3) gauge theory Jl9|] suggests that this shift 
corresponds to a decrease of the critical temperature of about 15%. 

The onset of deconfinement with increasing temperature is reflected in a sudden 
rise of bulk thermodynamic observables. For B > we thus expect to find a shift 
of this onset region to smaller temperatures. Otherwise, however, we expect that 
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Figure 4: Polyakov loop susceptibility (xl) and derivative of the Polyakov loop with 
respect to (3 (xp) for different values of B and lattices with spatial extent N a = 8 
and 10. For B = we only show the Ferrenberg-Swendsen interpolations of the data 
for better visibility. 
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Figure 5: Plaquette expectation value for different values of B and lattices with 
spatial extent N a = 8 and 10. The solid line shows a spline interpolation for the 
zero temperature plaquette expectation value, P , calculated on an 8 4 lattice. 
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observables like e.g. the free energy density show a temperature dependence similar 
to that in the pure gauge sector. In the high temperature limit we expect to find a 
gluon gas slightly modified due to the presence of static quarks. 

The free energy density can be calculated following the same approach used for 
B = 0, i.e. through an integration over differences of plaquette expectation values 
calculated on asymmetric (8 3 x 2) and symmetric (8 4 ) lattices jnj. In Figure |5| we 
show the plaquette expectation value, 

1 (S G ) , (3.12) 



6N%N T 

calculated for different values of B. The behaviour is similar to that of the Polyakov 
loop; with increasing ub the transition region broadens. The area between these 
data and the corresponding zero temperature results (solid line) increases. From 
this we obtain the free energy density^, 

/ p 



-6iV T 4 / d/3'[P -P] . (3.13) 

J0n 



Po Jp 



As can be seen from Figure |6| the free energy density decreases at fixed temperature 
with increasing B. This results in a shift of the onset of deconfinement to smaller 
temperatures. In the high temperature limit we indeed find that the free energy 
density is close to that of an ideal gluon gas at B — 0. The contribution of static 
quarks remains small. 



4 Conclusions 

We have formulated the quenched limit of QCD at non-zero baryon number. Al- 
though this formulation still leads to a path integral representation of the partition 
function with an integrand that is not strictly positive it can be handled quite well 
numerically for moderate values of the baryon number. A first numerical simulation 
of the gluon thermodynamics in the background of static quarks shows the expected 
behaviour. We find indications for a region of coexisting phases, which broadens 
with increasing baryon number density. The critical temperature, at which the 
system enters this region, is shifted towards smaller temperatures with increasing 
baryon number density. At high temperature we recover the physics of a gluon gas 
similar to the B = case. 

°We define as zero temperature lattice a symmetric lattice, iV r = N a . The plaquette expectation 
value (Po) calculated on the zero temperature lattice is used for normalization of the free energy. 
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Figure 6: The negative of the free energy density in units of T 4 evaluated for different 
values of B and lattices with spatial extent N a = 8. The abscissa is a linear scale in 
terms of the gauge coupling (3 covering the interval [4.74, 6.06]. We only show tick 
marks at the critical couplings on lattices with temporal extent iV T = 2, 3, 4, 6 and 
8, which correspond here to temperatures T/T c = 1, 1.5, 2, 3 and 4. 

Although the transition clearly reflects the physics of deconfmement it may not 
be a true thermal phase transition that leads to discontinuous changes in thermo- 
dynamic observables. A clarification of this point will require calculations on larger 
lattices. If the statement turns out to be correct one may have to think about the 
deconfmement transition at non-zero baryon number in terms of a non-thermal per- 
colation transition which anyhow seems to be a more appropriate physical picture 
for the deconfmement transition at zero temperature [20|j . 

We also find (indirect) evidence that the heavy quark potential does tend to a 
finite value at large distances already in the hadronic phase. The increase of the 
Polyakov loop expectation value with increasing baryon number density suggests 
that already at low temperatures string breaking starts at short distances. The 
heavy quark potential thus may be significantly modified in dense nuclear matter. 

The formulation of quenched QCD in the presence of static quarks seems to be 
an appropriate model for the thermodynamics of QCD at non-zero baryon number. 
It may be similarly useful for the analysis of thermal properties of hadronic matter 
at non-zero density as the pure SU (3) gauge theory has been for the understanding 
of the finite temperature deconfinement transition. 
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Temperature Phase Transitions in Particle Physics, EU contract no. ERBFMRX- 
CT97-0122 and the DFG through grant no. KA 1198/4-1. 
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A Appendix 



In this appendix we will derive the representation of fs given in Eq. 2.19 and give 
the explicit representation of fs for a few values of B. 

The starting point for deriving Eq. |2.19| is the representation of fs given in 
Eq. [HI, 

f b = (2K) BN ^f[D nh (x k ) , (A.l) 

X k=l 

where B' < B is the number of distinct sites x k appearing in X and n k is the quark 
occupation number for these distinct sites x k . They obey the constraint J2k n k = B 
with < n k < 6n/. 

Let us set up some general rules for the evaluation of fg. Any term in the sum 



appearing in Eq. |A.1| is characterized by Qnf numbers k = (hi, k 6nf ) where k{ 
indicates how many sites are occupied z-times. The numbers k{ are constrained 
by Hiiki = B. We thus may replace the vector X = (xi,...,xb) by a vector 
Y = (x*x, xb>) which only lists the distinct sites and is ordered according to the 
number of times these sites appear in X. The additional vector k keeps track of this 
information, i.e. the first k\ sites in Y appear only once in X, from k± + 1 to hi + k<i 
we label the sites which appear twice in X and so on. We then may write 



f B = (2K) BN ^S k (Y)^(2K) BN ^Sk ■ (A.2) 

Y,k k 



with 



6rc/ hi 

Sk(Y)=S (kl> ... Mnf) (Y) = l[ J] A&) • (A.3) 

i=l j=hi-i+l 

Here h = and hi = J2i<j<i % an d thus h &tlf = B' . Now we can perform the sum 
over Y for a fixed set of occupation numbers for the B' distinct sites, which is given 
by the vector k, 



S k = Y,'S k (Y) . (A.4) 

Y 
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The prime on the sum reminds us that in doing this sum we have to avoid, of course, 
double counting. A pair of sites Xi ^ Xj which have identical occupation numbers 
should appear only once in the sum, i.e. interchanging X{ and x 3 - should be counted 
as one configuration. This means that all Y which only differ by a permutation 
within an interval (hi + 1, h i+ i) are equivalent. We can give up this restriction and 
divide by appropriate factors, 

S k = (HkA . (A.5) 

Now we only have to insure that all Xi ^ Xj. We also want to eliminate this 
restriction and do independent sums over all Xi. This is achieved by performing 
one of the sums over Xi without any restriction and at the same time correct for the 
constraint by subtracting a term where two summation indices have been contracted. 
This process of factorization and contraction is repeated (B' — 1) times. Let us 
illustrate this by doing the first step explicitly: 



Y X\^X2 — ^S B l j=l 

LS' 

(X 



= (£Au(*i))(^£ 



3> 

xi 7 K X2-..^x B , j=2 

B f B' 

£ £ D ni (xt)T{D n] (x 3 ) 

t=2 x 2 ...^x B , 3=2 

= [A*] £ f[^n,(^)-£ £ D ni (x t )\{D n3 (x 3 ) , 

X2-..^x B / j=2 t=2 x2...^x B , j=2 

where [...] denotes the sum over the lattice taken on a single configuration, i.e. 

I1AJ =EIIAh(*) • ( A - 6 ) 

i xi 

When continuing this process of factorization and contraction we arrive at a cluster 
decomposition of the B' factors we started with. Contracting two summation indices 
leads to a factor (—1). We thus obtain a representation of as a sum over products 
of clusters 



r 6n f n«i 

F(a) = F(a u a 6n/ ) = (-l)^(j — 1)' II V 

L i=i 



(A.7) 
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with 



6rif 

j = E°i • ( A - 8 ) 
1=1 

Here we have explicitly given the combinatorial factor for a cluster of length j, i.e. a 
factor (—1) for each of the (J — l) contractions needed to generate a cluster of length 
j and a factor (j — l)! for the number of ways one can create such a cluster. We also 
included a combinatorial factor that takes into account that the permutation of ai 
factors Di appearing in a cluster does not lead to a new cluster decomposition. 

In total there are k[\ possibilities to distribute the various terms D[. However, 
we have to take into account that each cluster can appear several times, i.e. its 
degeneracy is g( ai ,..., o 6n )■ Also the permutation of identical clusters does not lead 
to a new cluster decomposition. We thus arrive at the representation 

T,s k (Y) = (YikAY,s({g a })U(—^) 9a ) ■ (A.9) 

Y \=l J {g a } {a} K 9a- t 

where a = (aj, ... a§ nf ) and the product runs over the set of vectors a which satisfy 

6rif 

£^<s ■ ( A - 10 ) 



Associated with each contributing vector a is a sum over the non-negative integers 
g a which is symbolically represented in ( |A.9| ) by the sum over {g a }, 



B int(S/2) i(a) 

£ = E £...£.-, (A.n) 

{g a } 9(i,o,. ..,o)=o 3(2,. ..,o)=o g a =0 

where i(a) denotes the integer part of B/Ya=i ^ a i- The Kronecker-5, 5({g a }), ap- 
pearing in ( |A.9| ) summarizes the constraints the set of summation indices has to 
satisfy, 

6n/ 6n / 

5({U)=^ [B-J2 lk i )Ti 6 (ki-Y,^) . (A.12) 

We note that the combinatorial factor appearing in (|A.9|) in front of the sums will 
just cancel the factor appearing in (|A.5 ); in fact, they do no longer depend on the 
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k 


S k {Y) 




001000 




1 


110000 


D 1 (xi)D 2 (x 2 ) 


1 


300000 


D 1 (x l )D 1 (x 2 )D 1 {x 3 ) 


3! 



Table 1: 3-quark contributions 

actual choice of the vector k. We thus may easily sum over all possible vectors k 
which yields essentially the result given in (|A.9| ) with a less restrictive constraint on 
the possible choice of the summation indices {g a }, 

/ B ,(2Kr^i(5-^iJn(^(«H , (A.i3) 

{g a } V 1=1 {a} J {a} K 9a- J 

where the allowed set of numbers {a} is constrained only by ( |A.10|) . 

This is the representation of Jb given in Eq. [2.19| . In the following we will give 
the explicit representation for a few small values of B. 

A.l B=3 

We get 3 contributions to f 3 : 

/ 3 = (2k) 3Nt ^Sooiooo + Snoooo + •S'sooooo^ • (A. 14) 

The terms are sums over products of terms Di which are defined in Eq. |2.14 . 
The three different contributions are given in Table [p. We thus obtain 

A = (2^) 3 ^([^3] - [DA] + Im + [Di\{[D 2 ] - \[Dl\) + ifA] 3 ) .(A.15) 
A.2 B=6 

For B = 6 one gets 11 contributions to f%: 



k 
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(A.16) 



k 


S k {Y) 




000001 


D 6 (xi) 


l 


100010 


D l (x 1 )D 5 (x 2 ) 


l 


010100 


D 2 {xi)D A (x 2 ) 


l 


002000 


D 3 (x 1 )D 3 (x 2 ) 


2! 


200010 


Dx{x 1 )D 1 {x 2 )D i {x 3 ) 


2! 


111000 


D 1 {x 1 )D 2 {x 2 )D 3 {x 3 ) 


1 


030000 


D 2 (xi)D 2 (x 2 )D 2 (x 3 ) 


3! 


301000 


D 1 {x 1 )D 1 {x 2 )D x {x 3 )D 3 {x^) 


3! 


220000 


D 1 {x 1 )D 1 {x 2 )D 2 {x 3 )D 2 {x A ) 


2! 2! 


410000 


Di(xi)D 1 (x 2 )D 1 (x 3 )D 1 (x 4 )D 2 (x 5 ) 


4! 


600000 


D 1 {x 1 )D 1 (x 2 )D 1 (x 3 )D 1 (x i )D 1 (x 5 )D 1 (x 6 ) 


6! 



Table 2: 6-quark contributions 

The different 6- vectors k and the corresponding contributions S).(Y) are listed in 
the Table |. 

This again gives rise to contributions, which can be ordered according to the 
number of clusters contributing 

; 6 = (2k) 6 ^]T> 6 (z) , (A.17) 

i=l 

with 

« 6 (1) = —[Df] + [DiD 2 ] - *-[DlDl\ + l -[Dl\ - [D\D 3 \ + 

2[ J D lJ D 2J D 3 ] - l -[Dl\ + [D*D 4 ] - [D 2 D 4 ] - [D 1 D 5 ] + [As] (A.18) 

ae(2) = -^mm] - 6[AA.]) + \[D,D 2 f + 
[m\[Dl\-[DiD 2 ]+\[Dz])- 

1 -{2[D 2 ] - [Dl]){\[Dl] + \[D\] - [D\D 2 ] + [AA,] - + 

[AK^^f] - [D\D 2 ] + [A^ 2 ] + [£?Ai] - [D 2 D 3 ] - 
[A A] + [A]) (A.19) 
« 6 (3) = ^(2[A]-[^ 2 ]) 3 + 
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\[Di\ 2 {\[Dl\ + \[Di] - [D\D 2 ] + [DA\ - [-D4]) (A.20) 
«e(4) = ^i] 2 (2[£ 2 ] - [D 2 ]) 2 + |[^i] s (|[^?] - [DA] + [A*]) (A.21) 
a 8 (5) = ^[AlWd - [JJ 2 ]) (A.22) 
06(6) = ^[A] 6 • (A.23) 



These coefficients have been generated with Mathematica. In total there are 58 
different terms contributing to fa. It is apparent that the number of terms increases 
rapidly with B. For B — 12 and rif — 1 there are 58 different vectors k which give 
rise to 2739 terms in f 12 . 
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